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Abstract. We study synchronization dynamics in networks of coupled oscillators with bimodal distribution 
of natural frequencies. This setup can be interpreted as a simple model of frequency synchronization dy- 
namics among generators and loads working in a power network. We derive the minimum coupling strength 
required to ensure global frequency synchronization. This threshold value can be efficiently found by solv- 
' ing a binary optimization problem, even for large networks. In order to validate our procedure, we compare 

its results with numerical simulations on a realistic network describing the European interconnected high- 
CZ2 . voltage electricity system, finding a very good agreement. Our synchronization threshold can be used to 

r/3 1 test the stability of frequency synchronization to link removals. As the threshold value changes only in 

O , very few cases when aplied to the European realistic network, we conclude that network is resilient in this 

" 55 ■ regard. Since the threshold calculation depends on the local connectivity, it can also be used to identify 

critical network partitions acting as synchronization bottlenecks. In our stability experiments we observe 
<""j , that when a link removal triggers a change in the critical partition, its limits tend to converge to national 

Ql* borders. This phenomenon, which can have important consequences to synchronization dynamics in case 

of cascading failure, signals the influence of the uncomplete topological integration of national power grids 
at the European scale. 

' PACS. 05.45.Xt Synchronization; coupled oscillators - 89.75.Fb Structures and organization in complex 

systems 

p: 

|- |. 1 Introduction Power systems are formed by a large number of genera- 

" tors interconnected in a complex pattern to supply energy 

(N The electrical power grid is an example of complex system to final consumers. Modeling the complete set of variables 
relying on the proper interaction between a great num- tnat characterizes the whole system represents a tanta- 
ber of elements [J. Recently it attained considerable at- lizin S effort - Here we follow the direction to reduce the 
JTj ] tention of the complex systems community due to many complexity of individual elements in favor of the complcx- 
rS ', features as, for instance, its non-trivial topological prop- ^ °f tne interaction pattern. 

^ ■ erties [21314) . the presence of cascading effects [516171819] . Our point tries to be a compromise of two opposite 
self-organized criticality as a possible explanation for the directions of research. On the one hand a more technical 
frequency of blackouts [TU] , its interaction with other net- perspective (electrical engineering) relies on the compli- 
work systems |11I12| . and also synchronization phenom- catcd (not complex) arrangement of primary, secondary, 
ena |13ll4j . Besides, current organizational trends are pos- and tertiary circuits to analyze the frequency stability of 
ing new research challenges related to the control of elec- the generators |17ll8j . On the other hand, a more statis- 
tic power grids. Two main prevaliling tendencies are the tical view (network science) analyzes the global behavior 
progressive integration of national networks into a European- mainly from the topological static features of the network; 
wide one [T5], and the shift from centralized energy pro- hence, the role of the different elements that form the net- 
duction towards more decentralized smart grids [TB]. works (nodes for generators and loads, links for distribu- 

tion lines) is uniquely associated with topological mea- 

a e-mail: slozanoOiphes . cat sures as, for instance, different types of centrality (local 

b e-mail: buzna@frdsa.uniza.sk or global) |19I20| . As discussed in [3T] these two lines of 

c e-mail: albert.diaz@ub.edu research can usefully inform each other in order to bring 
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new progress into an integrated study of very large sys- 
tems formed by units with internal complex dynamics. 

This paper can be also seen as an extension of the 
work presented in |14| , since here we are providing a gen- 
eral formalization of the synchronization threshold where 
the previous relation was a special case. Moreover, in the 
present work wc apply this generalized approach to study 
the influence of topological features on the stability of fre- 
quency synchronization in a realistic system. 

In section [5] we start the paper by reviewing the dy- 
namical description of power systems as networks of Ku- 
ramoto oscillators [35]. In section [3] wc derive a new for- 
mula describing the critical coupling strength and test its 
validity using realistic network topologies. In section 3] we 
analyse the sensitivity of frequency synchronization to link 
removals, and our findings and conclusions are summa- 
rized in section [SJ 



2 From swing and flow equations to 
oscillatory dynamics 

Our aim is precisely to incorporate simple dynamical be- 
havior of generators, i.e. not to deal with a complete dy- 
namical description involving electromagnetic fields and 
flow equations, but considering only the basic mechanism 
that ensures one of the crucial dynamical properties of the 
power distribution system as a whole, namely its synchro- 
nization. Synchronization is understood as the ability of 
an (unsupervised) system to keep a global state where gen- 
erators and loads (with intrinsic different frequencies) run 
with the same effective frequency (that of the distribution 
system, 50/60 Hz). To this end, generators and loads are 
described as phase oscillators in terms of "swing" equa- 
tions. As it has been shown recently, this equation can be 
derived from a proper energy balance in the generator |13j . 
The mechanical/thermal power brought to the generator 
gives rise to different contributions 

Psource — Pdissipated ~t" P accumulated ~t~ Ptransmitted (1) 



where 



Pdissipated 



(2) 



corresponds to the rate at which energy is dissipated due 
to the rotation of the mechanical rotors (turbines) and 7 
is a damping coefficient, 



r accumulated — 1 7 . " — luu 

at 



(3) 



corresponds to the rate at which kinetic energy is accumu- 
lated, with I being the inertia moment of the mechanical 
rotor. In Eqs. ©-(0 6 is the angle of the mechanical ro- 
tor. Finally, the energy is transmitted to the distribution 
system through the lines. The transmitted power is pro- 
portional to the sinus of the phase difference between the 
voltages of the elements at the two end-points of the line. 
Here it lies one of the crucial points of the mechanical- 
electrical approach, i.e. the assumption that the voltage 



angles and the rotor angles are the same. Thus, we can 
write the power transmitted from element i to j 



P 



-R 



MAX , 



(4) 



with P^f AX being the maximum power transferred be- 



tween generators i and j 
each element 



Hence, we can write for 



P 



1$ 



dt 
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R 



MAX , 



i) (5) 



assuming that Plf AX = if elements i and j are not phys- 
ically connected. This equation can describe generators as 
well as loads with the only difference that for loads Pi < 0. 
The frequencies of the different elements cannot be very 
far from the standard frequency of the distribution system 
(i?=50/60 Hz), then we write 



= Qt 



(6) 



assuming (pi <C Q. Inserting this expression in ((5]) and 
keeping linear terms in tf>i we get 



Pi 
2jiO 



n 

~2 



A 



■'ri 



which, for simplicity, can be recast as 



x S m(i Pj -i Pi ) 
(T) 



E 1 

3 



sin(y>j -ifi), 



(8) 



having introduced uii as the natural frequency of unit i and 
Wij the coupling term. The simple form of Eq. (|8) recalls 
the swing equation and this is the reason for its name in 
the technical literature, but the forcing term is what links 
the evolution of the phases of the element to the distribu- 
tion system. In terms of dynamical systems and complex 
networks the meaning of Eq. ([8J is obvious. It represents 
the evolution of phase oscillators with an additional inertia 
term |24I25| embedded in a complex network whose 
elements interact in a weighted way through the sinus of 
the phase differences between the end points of a physical 
distribution line. 

The Kuramoto model with inertia has been an- 
alyzed previously in the literature and its main 
effect is on the type of transition that changes 
from second to first order when the coupling is 
increased. Although in |24| it is reported that in- 
ertia modifies the phase diagram of the system, 
one has to bear in mind that in that case the 
authors consider also the effects of noise in the 
system. For noiseless systems, as those considered 
in 1 27 ,28 1, it is shown that there are two critical 



values of the coupling strength and a hysteresis 
behavior within this range of values. The lower 
value, below which the incoherent state prevails, 
is not modified with respect to the original Ku- 
ramoto model; but, the upper value, above which 
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the synchronized state is the only possible solu- 
tion, changes with the inertia contribution. In the 
intermediate region, there is a coexistence of the 
synchronized and incoherent states. The hystere- 
sis behavior shows that the stability of the syn- 
chronized state is broken at the lower value of the 
coupling strength, and this is the important ef- 
fect we will look at in this work. Under usual op- 
erational conditions the system is frequency (not 
phase) synchronized and we want to analyze how 
any disturbance (removal of some links) can affect 
this state, this means that only the lower value 
of the coupling strength matters, since it is the 
value below which the synchronized state is not 
achievable. Additionally, it has been shown [29 30j 
that the inertia term modifies the time scale for 
the achievement of the synchronized state but not 
the existence of the state. For all these reason we 
will remove the inertia term, which only speeds 
up convergence of the numerical simulations but 
it does not affect the conditions under which the 
system can achieve synchronization. 



3 Stability of frequency synchronization 

What makes a systems of oscillators interesting is its dy- 
namical behavior and, from a collective point of view, the 
most relevant feature is its ability to synchronize. As we 
will show in this section, for a population of non-identical 
oscillators, perfect phase synchronization is not possible, 
whereas frequency synchronization requires a minimum 
coupling strength. 

3.1 Estimation of the minimum coupling required 

As it was shown in [14] , when we neglect the inertia term 
from Eq. (|8|). we can relate the synchronization threshold 
to the network topology. We consider the networked pop- 
ulation of oscillating units, organized in a graph Q com- 
posed from a set of nodes N and a set of links L. Each 
unit (node) i G N is characterized by a natural frequency 
u>i and a phase angle ipi. The dynamics of these units is 
then governed by 

ipi = ojj + a ay sin(</?j - cpi), (9) 

3 

where ay is an element of the adjacency matrix which 
takes value 1 when nodes i and j are connected by a link 
and value otherwise. Note that compared to Eq. we 
have neglected the inertia term, we assume all links to be 
identical, and their coupling strength is characterized by 
a parameter a. Motivated by the existence of generators 
and loads in power systems, we split the nodes into two 
populations of identical elements, where the members of 
the first population N + C N behave like power produc- 
ing units, and thus w,; = w+ > for i G N + . And vice 
versa, the members of the second population A_ C N, 



where A_ U N+ = N and iV_ D N + = 0, behave as power 
consuming units and Ui = uj- < for i G JV_. Since the 
number of generators and the number of loads are not 
necessarily equal, we normalize u>i values in the following 
way u> + = 1/|A + | and w_ = — l/|iV_|. The choice of the 
positive sign for the natural frequency of the generators is 
not arbitrary; note that in eq. (O this is the case whereas 
for loads this term is clearly negative. 

In general, for a population of Kuramoto oscillators 
there is a balance between the two terms in eq. ©• On the 
one hand, if there is no coupling all units follow their nat- 
ural frequencies a;, and there is no frequency synchroniza- 
tion. On the other hand, when the coupling is very large 
the effective frequencies ipi tend to zero. Then, there will 
be a critical value of the coupling strength a above which 
the system synchronizes in the sense that the effective fre- 
quencies become equal. Note that this frequency synchro- 
nization docs not imply phase synchronization, which is 
only possible, for a population of non-identical oscillators, 
when the coupling strength goes to infinite. 

Following the derivations presented in |14j . we can 
write for any unit in the synchronized state 

= uji + a a,ij sin(<y9j — ipi). (10) 

Here it is easy to see that a necessary, but not sufficient, 
condition for the natural frequency of unit i being is 
that 

'>°t=^ (11) 

Ki 

where ki is the degree of the node i. From this expression 
it is, in principle, possible to find a global bound for the 
coupling strength. This bound can be, however, far from 
the real value, since there are some geometric constraints 
[31j (phase differences that cannot be maximized simulta- 
neously) involving more than one unit. Actually, we can 
write as many equations as partitions of the network into 
two non-overlapping clusters, just by summing the equa- 
tions for the nodes and using the fact that the interaction 
is an odd function. Identifying the network partition as a 
set of nodes S, we can write for the sum of the Eqs. of the 
nodes in S 

= ^2 U i + ° X! a « ^(fj ~ Vi)- (12) 
i€S ieS.j^S 

from which we can find the critical values of the coupling 
for the partition S 



which generalizes (|lip , containing it as particular cases of 
clusters formed by individual nodes. Note that in [T5] we 
restricted ourselves to clusters of the same type of node, 
which gives good estimations for homogeneous distribu- 
tions, but not for a more even distribution as is the current 
case of interest. The goal is to obtain the global minimum 
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Fig. 1. Topology of the approximate model of European 
interconnected power system [9]. White circles represent 
positions of active (power injecting) stations and blue dots 
stand for passive (power consuming) stations. 



that ensures the existence of a synchronized state 



a > max 



E 



(14) 



3.2 Testing numerically the critical coupling 

To test our hypothesis we use the approximative data de- 
scribing the topology of the European power transmission 
network [5] composed from sixteen national subnetworks. 
The network topology was derived from maps published 
by UCTE [32] ■ Functional parameters, as for example po- 
sition of power producing and power consuming units, 
their capacities and capacities of power lines, were esti- 
mated in order to fit the publicly known volumes of cross 
border flows as precisely as possible. 

The network is depicted in Fig. [TJ where white circles 
represent positions of active (power injecting) units and 
blue dots stand for passive (power consuming) units (can 
be seen online when the figure is sufficiently zoomed) . The 
first two columns in Table [1] summarize the number of ver- 
tices and edges for the European network as a whole, and 
for each one of the national subnetworks (i.e. those parts 
of the network contained within each country boundaries) . 

According to the previous paragraphs, we have applied 
Eq. (fT4"|) to each one of the networks in order to deter- 
mine their critical coupling value a c (i.e. the minimum one 



needed to have frequency synchronization Jj. Specifically, 
we have applied a semi-analytical approach to determine 
the partition within the network maximizing the expres- 
sion in Eq. (|14j) . This procedure provides lower and upper 
bounds for a c , a l A and a A , respectively. The third and 
forth columns in Table [T] present a l A and <r A for each one 
of the considered networks. More details on the method- 
ology can be found in the Appendix. 

To verify the goodness of these values as approxima- 
tions to the critical coupling, we have to check whether 
they arc the minimum value of a allowing full synchroniza- 
tion of the system. We have done so numerically by sim- 
ulating frequency synchronization dynamics. More con- 
cretely, we have run synchronization dynamics described 
in Eq. @ and measured frequency dispersion (?*), which 
is the order parameter proposed in P3] to measure of 
the effective frequency dispersion: 



(15) 



ieN 



Fig. [U presents time evolution of the order parameter 
r for different values of the coupling a for the European 
network. For values far enough from the critical one, the 
system cither fluctuates steadily around a certain value 
(see a = 0.005 in the figure) or relaxes towards full fre- 
quency synchronization (a = 0.1 case). As we get closer to 
the critical value, we observe an initial tendency toward 
synchronization that is sharply broken by very strong fluc- 
tuations (er = 0.011). 

For each network we have executed an iterative pro- 
cedure running the dynamics with several coupling val- 
ues and checking whether r was relaxing towards 0. In 
each case, the minimum coupling satisfying this condition 
, , has been determined with a precision of four decimal 
places. 

Finally, a 1 a values in Table Q] (obtained analyt- 
ically applying Eq. ([14])) can be checked by com- 
paring them with these o^'s obtained numerically. 
Fig. [3] presents all values together to allow for a vi- 
sual comparison. Generally speaking, we observe a re- 
markable agreement between the two sets of results. Even 
for the three cases where the difference among analytical 
and numerical approximations is visible (i.e. Europe as a 
whole, Germany and Poland), the obtained values are of 
the same magnitude. 

This outcome supports our claim that Eq. (|14j) is a 
good estimation of the critical coupling c c , and opens the 
door to simple studies on frequency synchronization in 
networks. In particular, Eq. ([T4"]) can be applied to the 
study of frequency synchronization's stability against per- 
turbations. 



1 Absolute values of a c are fully dependent on u)+ and 
The values of lj+ and cj_ were chosen in a way that they sum up 
to 1 and -1, respectively. As the number of nodes in national 
networks differs it makes the comparison between countries 
difficult. 
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a = 0.005 
= 0.011 
a = 0.012 

= 0.1 




le+05 



5e+04 

Time-steps 

Fig. 2. Temporal evolution of frequency synchronization 
in the European network for different coupling values. 
Frequency dispersion r (the order parameter introduced 
in [Tl]) is used to show the effect of coupling a below 
(0.005 and 0.011) and above (0.012 and 0.1) the critical 
value. Notice that the plot corresponding to a = 0.011 
presents a characteristic behavior of a system very close 
to criticality. First it relaxes towards synchronization but, 
at a certain point, it experiences sharp fluctuations that 
are stronger than those in the case a = 0.005. Simulations 
start with all phases set to so, as the dynamics are purely 
deterministic, only one realization per coupling value was 
needed. 
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Fig. 3. Comparison of a l A values determined by means of 
Eq. (j!4[) and erjv obtained numerically using the iterative 
approximation method. As a guide to the eye, we have 
added a dashed line corresponding to erjv = 0^4 • Clearly, 
there is a remarkable agreement. This result supports our 
claim that Eq. (|14[) is a good estimation of the minimum 
coupling needed to assure frequency synchronization. 



Country Name 


tvt 
In 


T 

L 




°A 


Europe 


1254 


1943 


0.00818157 


0.00818253 


Austria 


36 


429 


0.112499 


0.1125 


Belgium 


22 


21 


0.116666 


0.116667 


Croatia 


17 


20 


0.25 


0.250001 


Czech Republic 


34 


52 


0.090909 


0.09091 


Denmark 


8 


8 


0.466666 


0.466667 


France 


318 


519 


r\ r\ A a no 

0.034482 


r\ (~\ a /too 

0.034483 


Germany 


229 


313 


O.Uzy411o 


n noo /iioo 
U. 0294120 


Hungary 


27 


36 


0.166666 


0.166667 


iiaiy 




90/1 


U.U4: ( OU ( 4 




Luxemburg 


3 


2 


0.999999 


1 


Netherlands 


22 


24 


0.290598 


0.290599 


Poland 


99 


140 


0.043478 


0.043479 


Portugal 


24 


44 


0.188889 


0.18889 


Slovakia 


25 


30 


0.166666 


0.166667 


Slovenia 


8 


8 


0.333333 


0.333334 


Spain 


193 


316 


0.03125 


0.031251 


Switzerland 


47 


76 


0.115384 


0.115385 



Table 1. Network size and the values of a l A and a\ for the 
whole European high-voltage electrical network and the 
national networks. Values a l A and <j\ have been calculated 
according to Eq. (fl~4"| following the procedure described in 
the Appendix. 



4 Resilience of power systems in terms of 
frequency synchronization stability 

There is a bulk of literature analyzing the robustness of 
power systems. Since power systems can be represented as 
networks, most of these approaches deal with the robust- 
ness of the system to the removal of links or vertices. This 
kind of robustness analysis has been addressed in many 
different ways. Initial purely topological approaches |33j 
have shown to be limited and recent ones combine 
structural changes with different processes such as flow 
dynamics |6l34j . 

Following up from this literature, we address the ro- 
bustness of power systems to link removal from the view- 
point of the stability of frequency synchronization. A straight- 
forward way to make such an analysis is to calculate to 
what extent the removal of a link from the network mod- 
ifies the critical value of the coupling a c for the whole 
system. 

More concretely, the idea is to test whether remov- 
ing certain links from the network can increase the min- 
imum coupling value required to ensure synchronization 
(i.e. making it less stable) or the other way around, de- 
crease it (subsequently improving its robustness). 



4.1 Effect of link removal on synchronization stability 

Our starting point is the complete European network 
Fig.Ushows the location of the partition maximizing Eq. (TT 
and, therefore, determining a c . Then, for each one of the 
links in the European network e, we take it out from the 
grid, calculate the new value of a l A (e), and restore it to its 
position. 
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Fig. 4. Critical partition of the European network. The 
shadowed area corresponds to the partition defining the 
critical value of coupling (a l A ) when the European net- 
work is complete, and the links marked in green connect 
it to the rest of the network. Notice that, given the condi- 
tions stated in section 2, the complementary partition (i.e. 
the whole network except the shadowed area) presents the 
same critical coupling, so any of the two could be consid- 
ered the critical one. 
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Fig. 5. Impact over global frequency synchronizity of re- 
moving each edge of the European Network. New critical 
coupling values a 1 a were obtained from Eq. ([14]) as de- 
scribed in the Appendix, and edges are ordered according 
to their Edge Bctwccnness Cs(e). Inset: Cumulative prob- 
ability distribution of a 1 a values. The network is robust 
to removal of all but very few edges, whose actual impact 
when taken out is not correlated with their centrality. 



As shown in Fig. [5J only few links in the whole Euro- 
pean network have an impact on the global a c value when 
removed. Among them, a minority decrease <j\ (i.e. make 
overall frequency synchronization easier if erased), while 
the other ones exercise the opposite influence, since their 
removal increases a l A value. 

We can also observe that such effects of link deletion 
are not related to link betweenness [35] . This result is es- 
pecially relevant, since purely topological measures (such 
as link betweenness or degree centrality) have been tradi- 
tionally used to assess resilience in power grids and other 
networked systems. On the contrary, here a strictly topo- 
logical approach is no longer valid and we need to include 
dynamical aspects in our analysis. Accordingly, under- 
standing the effects shown in Fig. [5] requires interpreting 
how a link removal affects synchronization dynamics in 
Eq. © and, as a consequence, modifies the global critical 
coupling. 

Fig. HI provides such a detailed view of the local struc- 
ture for five of the links which deletion changes the global 
critical coupling. Removed links are marked in red, the 
shadowed areas correspond to the partition becoming crit- 
ical after the removal, and green ties are their remaining 
connections with the rest of the network. 

Figs.EI a )-(d) show cases where the deletion of the link 
marked in red makes the critical coupling er c to increase. 
It is quite intuitive to see that the removed links were 
partition's 'boundary' links (i.e. pointing to outside of the 
partitions) and, therefore, their removal reduces the con- 
nections of the partitions to the rest of the network (e.g. in 
casejnja), for instance, from two to just the one marked in 
green). In accordance with Eq. (|14j) . this reduction leads 
to an increase of the partition's critical coupling, which 
eventually becomes the new global a c . 

Fig. Etc) presents an example of the opposite case (i.e. 
a c diminishing). Here the erased link was an inner one, 
connecting a leaf (i.e. a node of degree one) to the rest of 
the partition. Therefore its deletion reduces the size of the 
partition and, again according to Eq. (|14[) , makes the crit- 
ical value of the coupling for the partition smaller. This 
reduction then let other partitions in the network with a 
critical coupling lower than the previous global a c become 
the new critical partition. In this particular case, the new 
critical partition corresponds to the Iberian Peninsula. No- 
tice that in contrast with case Eld) , which presents the 
same resulting critical partition, in this case the number 
of remaining boundary links are 3, so the critical coupling 
of the partition is lower. 



5 Conclusions 

We studied how synchronization dynamics of networked 
oscillators depend on their topological configurations. When 
neglecting control mechanisms present in real power sys- 
tems and considering only some aspects of the swing and 
flow equations, oscillatory dynamics can be captured by 
the Kuramoto model assuming a bipolar distribution of 
natural frequencies. For this setup we are here giving accu- 
rate analytical estimation of the critical coupling strength. 
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Fig. 6. Examples of link removals modifying the critical 
value of coupling a\ of the European network. In cases 
a) - d), the link removal increases a l A , while it is reduced 
in case e). For all cases, removed links are marked in red 
and indicated with an arrow, the new critical partition is 
marked as a shadowed area and the links connecting it to 
the rest of the network are marked in green. More details 
are provided in the main text. 



It leads to the solving of a combinatorial optimization 
problem on the graph. This can be efficiently done us- 
ing methods of binary linear programming even for very 
large networks with an arbitrary precision. We validated 
our results on realistic network topologies and by compar- 
ing them with numerical simulations we found very good 
agreement. 

Furthermore, we studied the robustness of the synchro- 
nization process with respect to removals of single links. 
When links are removed, only a small percentage of them 
have an impact on the synchronization threshold. Thus 



network is in this respect very robust. While some links 
are increasing the synchronization threshold the others are 
decreasing it. The concrete effect depends on whether they 
are decreasing the size of the critical cluster or affecting its 
connectivity with the rest of the network. We also tested 
whether the positions of links influencing the synchroniza- 
tion threshold is correlated with their global topological 
properties as is sometimes assumed when assessing the 
resilience of real-world network topologies. Here we find 
no correlation. This can be explained when interpreting 
the Eq. (|I4[) which estimates the synchronization thresh- 
old. For a given subset of nodes the threshold value is 
higher the higher is the sum of natural frequencies for 
given subset of nodes and the less tightly are these nodes 
connected with the rest of the network. Thus the synchro- 
nization threshold depends on the local properties of the 
critical network partition. 

When analyzing the geographical positioning of criti- 
cal partitions, we often found them being identical to bor- 
ders between European countries. For example, for the 
original network without removed links the critical clus- 
ters are connected along the western French border (see 
Fig. [4j. When links are removed, in some cases critical 
clusters are connected along the Spanish-French border 
(see Figs.lBJi andHjk). Another example is the critical par- 
tition corresponding to the Dutch network (see Fig. [6}d). 
These observations indicate the tendency of national net- 
works to be more tightly connected internally than across 
borders, which is a signature of an on-going (incomplete) 
European-wide integration process. This fact can have 
positive consequences from a practical point of view, as it 
leads to a natural tendency in the network to split along 
these areas in critical situations. The major disruption in 
the synchronization of the European network in Novem- 
ber 2006, which resulted in a temporal division of the sys- 
tem into three regions mainly following existing or former 
national boundaries, is an illustrative example of this de- 
pendence [56] , 

This paper is contributing to a better understanding of 
the interplay between the network topology determining 
the spatial positioning of network elements and frequency 
synchronization dynamics. Future steps could include the 
study of the cascading behaviours, whereas the Kuramoto 
equations allow to define link flows, or investigation of the 
validity of our results when considering a more realistic 
model, e.g. by introducing heterogeneous links or various 
types of generators and loads. Finally, this article was 
focused on the existence and stability of frequency 
synchronization, where the inertia has a limited ef- 
fect. Possible extensions of our work could address, 
for instance, the non-trivial influence of the inertia 
term on transient states pointed out in section 2. 
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Appendix: Methodology 



The problem of identifying the subset of nodes maximiz- 
ing Eq. (1141) can be efficiently solved by using binary linear 
programming |37j . First, we formulate it as a combinato- 
rial optimization problem: Let us have a graph Q (N, L) 
where N is a set of nodes and L is a set of links. Each 
node has a real value u>i associated to it. The problem is 
to find a nonempty subset of nodes S C N maximizing the 

expression ^^' es — — . 

Thus, we need to decide which nodes are, and which 
nodes are not, included in the subset S. Such decision can 
be modeled by a set of binary variables £ {0,1} for 
i £ N. Variable j/j = 1 iff i £ S and otherwise yi = 0. In 
the denominator of Eq. (|14p . we consider the number of 
links connecting nodes belonging to the set S with nodes 
which do not belong to the set S. To count them, we 
introduce for each link connecting nodes i and j a binary 
variable Xij £ {0, 1}. Variable Xij takes value 1 iff either 
i £ S and j £ N - S or j £ S and i £ N - S and 
Xi.j = otherwise. Then we can formally formulate this 
combinatorial optimization problem: 



Maximize / = 



Yl(i,j)eL X i,3 



subject to 



— Uj Hi 
^ Hi ~~ Uj 



(16) 



(17) 



for(i,j)£L (18) 
for £ L (19) 



Vi , x id £ {0, 1} for i £ N, £ L (20) 

Constraint (|17[) ensures that the set S is nonempty. 
The set of constraints (|18[) and (|19[) make sure that for 
each link connecting nodes i and j variable Xij = 1 if yi = 
1 and yj = or when j/i = and yj = 1. The objective 
function (fTB"]) is forcing variables Xi.j to take value of zero 
whenever it is possible. Therefore Xij is zero if yi = and 
Uj = or if yi = 1 and yj = 1. 

Note that we obtained an optimization problem with a 
non linear objective function (|16l) that is optimized with 
respect to the linear set of constraints (fT7|) - ([T9| . When 
substituting the value of the objective function ([TTJ)) by 
the variable a we can rewrite the problem ([T6l - (fT^)) as: 



Maximize f = o~ 
subject to 

(i,j)EL i£N 



(21) 
(22) 
(23) 



for £ L 
(24) 

Xi,j >V%- V] for (i,j) £ L 

(25) 

a > 0, Vi, Xi,j £ {0, 1} for i £ N, £ L 

(26) 

If we replace the variable a by an arbitrary constant 
c we obtain a linear optimization problem which can be 
easily solved by traditional integer solvers as, for example, 
XPRESS-IVE [38]. The remaining problem has a feasible 
solution only if c < a c and there is no solution if c > 
<t c . Thus, by using a binary search on c and repeatedly 
solving the optimization problem for different c values we 
can find an arbitrarily tight lower a l A and upper a\ bounds 
for a c . Moreover, values of variables yi corresponding to 
the lower bound a l A (when feasible solution exists) define 
which nodes belong to the partition whose a = a l A < a c < 
a A . 
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